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Numerical calculations of light propagation in random media demand the multiply scattered 
Stokes intensities to be written in a common fixed reference. A particularly useful way to perform 
automatically these basis transformations is to write the scattered intensities in the Chandrasekhar- 
Sekera representation. This representation produce side effects so that numerical tests are necessary 
to deal with the limiting situations of the small-particle (Rayleigh) and for/backwards scattering. 
Here a new set of basis functions is presented to describe the scattering of light by spherical par- 
ticles (Mie scattering) in the Chandrasekhar-Sekera representation. These basis functions can be 
implemented in a new algorithm to calculate the Mie scattering amplitudes, which leads straight- 
forwardly to all the scattering quantities. In contrast to the traditional implementation, this set 
of basis functions implies to natural numerical convergence to the above mentioned limiting cases, 
which are thoroughly discussed. 
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I. INTRODUCTION 

The scattering of an electromagnetic plane wave by 
a spherical (homogeneous, isotropic and optically linear 
material) particle of arbitrary size is known as Mie scat- 
tering QI H H Q . The numerical calculation of this scat- 
tering is relevant to several fields H H HH such as 
remote sensing (meteorological optics, radar detection of 
raindrops, lidar detection of clouds etc.), optical particle 
characterization (inverse problems) etc. More recently, 
magnetic Mie scattering has attracted the attention of 
researchers [II El El 111 El E| 

In optical dense media, the light scattered by a particle 
is sucessively rescattered and the electromagnetic wave 
phase coherence may not be totally destroyed by the par- 
ticle configuration avera ges producing interesting effects 
due to phase correlation ,16j . The role of numerical mul- 
tiple scattering description is twofold: it can be viewed as 
a tool in the comprehension of more fundamental aspects 
(memory of incident polarization state p^llMll^ . l20ll21| . 
weak localization pi I21 Ety) or as a tool for random 
media characterization (light scattering by b iological tis- 
sues [11 il il il |M 111 IM 111 S S 111, for in- 
stance). These procedures demand considerable numer- 
ical efforts and it is desirable to have very efficient and 
robust (covering a wide range of values in the parameter 
space) codes to perform such tasks. 

To deal with numerical multiple scattering simulation, 
basis transformations must be performed to represent the 
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local scattered Stokes intensities into a common fixed ba- 
sis, the laboratory basis. For instance, this procedure is 
successively repeated in multiple scattering of light in 
a Monte Carlo scheme. Other than being cumbersome, 
these basis transformations may increase the propagation 
of round-off errors in the numerical simulations and may 
ask for complementary tests in limiting situations. 

These basis transformations may be implicitly con- 
sidered writing all the wavevectors in the same labora- 
tory basis. This is the Chandrasekhar-Sekera represen- 
tation, which has been employed either in a radiative 
transfer equation calculation [3(| as well as in a Monte 
Carlo scheme 37]. The drawback of this representa- 
tion (not stressed in the literature) is the introduction 
of new difficulties to numerical calculation, notably the 
non-commutation of the two limiting cases: i) small par- 
ticle size compared to the wavelength (important to poly- 
dispersion calculations), and ii) for/backward scattering 
events (important to variance reduction in a Monte Carlo 
schemes and radar/lidar detection). 

In this paper the problem of limits non-commutability 
is pointed out and a new set of basis functions for Mie 
scattering calculation is proposed to the implementation 
of an algorithm. This algorithm naturally includes the 
small particle size and for/backward scattering. The pre- 
sentation is divided as follows. In Section [H] a brief re- 
view of Mie scattering in the scattering plane represen- 
tation is presented to set up the notation. Also, an orig- 
inal expansion of the scattering amplitudes up to fourth- 
order on the cosine of the scattering angle is calculated. 
In Section ITTT1 the Chandrasekhar-Sekera representation 
is reviewed. It is pointed out in Section IIVI the non- 
commutability of the limiting cases. In Section^ a new 
basis functions are presented to calculate Mie scattering 
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in the Chandrasekhar-Sekera representation and its im- 
plementation for numerical calculation is discussed. Con- 
cluding remarks are presented in Section IvTl 



II. MIE SCATTERING 

Let us consider a non-absorbing medium with (real) 
refractive index n m (n m 6 S), a sphere of radius a and 
(complex) refractive index n s (n s S C) to take emission 
or absorption into account. The origin of the labora- 
tory frame (i, y, z) is placed at the center of a spher- 
ical particle. An incident monochromatic (plane) wave, 
with electric field Eq and wavelength A (consequently the 
wavenumber k = 27r/A), propagates along the z-direction 

with wavevector fco = kz and it is scattered by the sphere. 
The interest is on the spherical scattered field E\ along 
the k direction, defined by the spherical angles 9 and <f> 
according to the laboratory frame. The scattering plane 
is formed by the vectors ko and k. Notice that the scatter- 
ing plane cannot be univocally defined for back/forward 
scattering cases and that the magnetic field (of the elec- 
tromagnetic wave) is not neglected, but it can always be 
obtained from the electric fields (because of the consid- 
eration of transverse wave). 



A. Scattering Plane Representation 

In a distance R from the center of the sphere, with 
(ka) 2 /kR <C 1 (far-field approximation), the scattered 
field is nearly transverse. The electric field lies on the 
plane orthogonal to k formed by the orthogonal directions 
8 and <fr, which form the spherical basis. Using i = \J — 1 
for the imaginary part, the scattered field in the spherical 
basis is written as: 



E„ 

E,h 



-ikR 



kR 



J(H,<p) 



E, r 

E,, 



(1) 



with /i = cos 6, 9 is the scattering angle and 4> is the 
azimuthal angle. The Jones matrix is the product of a 
diagonal matrix (because of the scatterer spherical sym- 
metry), known as scattering matrix with a rotation ma- 
trix: 



J(nA) 







S ± (jj) 



coscp 
— sin < 



sin ( 
cos ( 



(2) 



The rotation matrix projects the incident electric field 
(given in the laboratory frame) to the parallel and per- 
pendicular directions relative to the scattering plane. 
The scattering matrix then alters the field values via 
S»((jl) and S±(fi), which are the parallel and perpen- 
dicular scattering amplitudes (relative to the scattering 



plane), respectively, 

°° 2n + 1 

Swiv) = y^—f — rrrK^M + MnW], (3) 
' n(n + 1 

n— 1 v ' 

S±(fJ.) = — ~~rr [a>n'Kn (a*) + b n r n (/i)] , (4) 
■* — ' n(n + 1 

n— 1 v ' 



where a n and b„ 
and 

^n(p) 



are the Mic coefficients 1 
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with ir' n (n) — dTr n (p)/dn and P^ 1 ^) being the n tn Leg- 
endre polynomial of first order. In practice the summa- 
tions in n in Eqs. J2Jl and must be performed to 
Umax = fca + 4\/t^+2 jHUlj. 

The Mie scattering is dependent only on three quan- 
tities: the size parameter ka, (complex) relative refrac- 
tive index m — n s /n m (m € C) and rh = m/p,, where 
jl = [i s j[i rn is the relative ( s phe re/medium) complex 
magnetic permeability pH [U Hi El El d ■ The nu- 
merical calculation of this scattering event consists of two 
parts: one involving characteristics of the scatterer (Mie 
coefficients, which leads to several cross-sections) and the 
other involving the geometrical aspects (scattering am- 
plitudes, which leads to the phase function). 

The calculation of Mie coefficients a n and b n depend on 
ka, m and rh through the calculation of spherical Bessel, 
Neuman, and Hankcl functions and their first derivatives 
with respect to the argument (ka and mka). The diffi- 
culty in this calculation is that the recurrence relation- 
ships for complex arguments of the Bessel functions are 
not stable. This problem is solved writing the spherical 
functions in the ratio form [4(j for Mie coefficients and 
with the use of the continued fraction method developed 
by Lentz |4l) . 

If a desired precision in the results is known in ad- 
vance, it is possible to use the recurrences for a n and 
b n developed by Bohr en l42| and implement the scheme 
proposed by Cachorro Further improvements to this 
calculation have been compiled and studied in detail in 
Ref. 0. 

The second part of the calculation concerns the scat- 
tered field. It depends on the relative position of the 
detector with respect to the scatterer and source and 
also on the distance between scatterer and detector. The 
scattering amplitudes, which depend on the azimuthal 
angles, cosine of the scattering angle /1 and on the Mie 
coefficients, can be efficiently calculated using the algo- 
rithm created by Wiscombe |3^ | . 

Observe that the total scattering cross section is ob- 
tained throught the optical theorem, which is obtained 
expanding the scattering amplitudes up to second order 
on the scattering angle. 
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B. Scattering Amplitudes Expansion along the 
For/Backward Direction 

The expansions of the scattering amplitudes are pre- 
sented up to the fourth-order on the scattering angle. 



These expansions have been explicitly included since they 
are not easily found in the current literature. As shown in 
Ref. 0, p. 73], the Legendre polynomials can be written 
as: 



P„(cos6») = ^nCOS P 9 
p=S,S+2 

if n even 



5 = 



1 if n odd 



Mn = (-1)(-^ 2 



l-3-5---(n + p-l) 
2(»-p)/2[(n-p)/2]!p! ' 



(7) 

(8) 
(9) 



so that: Pn 1 \cos9) = sin 9dP n (cos 9)/d cos 9 and Eqs. JSJ and © are rewritten as: 7r„(cos#) = J2p=s s+2 t^nP coaP 1 ^ 
r n (cos#) = cos 6*7r„ (cos 0) — sin 2 9Tr' n (cos9) = Y^,p=s s+2 ^nP C0&p 2 #(1— psin 2 9). Expanding the functions up to 0(9 4 ) 
around the back and forward scattering directions, one obtains: 



7r„(oos0 - ±1) - + J H (±) _ K i±) 

r„(cos6»^±l) = ±H^± 



2 (±irv^ 
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No closed analytic form have been obtained for the latter summation. The scattering amplitudes, expanded up to 
0(9 4 ) around the back and forward scattering directions, are: 



S||(cos0-> ±1) 
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5(±1) 
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(14) 
(15) 

(16) 
(17) 
(18) 



From Eqs. (|14|) and l|15|) . one obtains the following 
properties for the exact forward and backward scatter- 
ing amplitudes: S±(l) = S\\{1), S±(-l) = -5||(-1). 
The limiting case cos 6* — s- 1 is used to calculate the total 
scattering cross section using the optical theorem while 
when cos# — > — 1 gives the backscattering (radar) cross 
section 0,HHI|. 



III. CHANDRASEKHAR-SEKERA 
REPRESENTATION 

With the scattering plane representation one can to 
express directly the photon "history" from the source 
to the detector. Nevertheless, for multiple scattering it 
presents some inconveniences. From one scattering event 
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to another, the electric field needs to be written in lo- 
cal basis (6, cf>, k), which depends on the position of the 
scatterer. For instance, to write the scattered field along 
a scattering sequence in a Monte Carlo scheme, the az- 
imuthal and scattering angles and all the distances along 
the sequence need to be recorded to perform the right ba- 
sis transformation back to the laboratory (fixed) frame. 
This bookkeeping is necessary when the reversed scat- 
tering sequences must be obtained. This representation 
becomes specially cumbersome for the calculation of the 
coherent backscattering enhancement (weak localization 

of light) El El El m. 

The difficulties presented above can be solved using the 
Chandrasekhar-Sekera representation 0, 0| . Consider 
th 

for example the k scattering event in a given scatter- 
ing sequence. As previously, an incident monochromatic 
(plane) wave with wavevector k K is considered and one is 
interested in the scattered field along the direction k K+ x- 
Note that in this case the incident direction is arbitrary, 
not necessarily along the z-direction, as in the scatter- 



ing plane representation. In a distance R K +x from the 
center of the sphere, with (ka) 2 /kR K+ i -C 1 (far-field 
approximation), the scattered field is nearly transverse 
and the electric field lies on the plane formed by the 
spherical basis vectors K +i and <t K .+x- The spherical 
basis (0 K , $ K , fc K ) is defined by the angle 9 K , which is 
obtained from the projection of the wavevector k K along 
the z-direction (cos K — k K ■ z) and by the azimuthal 
angle $ K , which is counted from x to the projection of 
k K onto xy— plane, relative to the laboratory frame. 

Equivalently to the scattering plane representation, 
one can write the scattered electric field E K +x along the 
directions (0 K +i, 3> K +i) as a function of the incident field 
E K along the directions (8 K ,$ K ): 

- g-ife-R K +i „ _ 
E K+1 (k K+ i) = -i — — J(k K+1 ,k K ) E K (k K ) , (19) 

where the Jones matrix is: 



J(k K+ x,k K ) = 
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(l,r) = 

(r,l) = 

(r,r) = 
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Xx ' 


"(J,r) 
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(r, r) (r, I) 




x 2 



(1 - - A+l) + fJ-nf^K+l cos(A$ K , K+ i) , 

-fi K sin(A$ KjK+ i) , 
fi K+ x sin(A<& K , K+ i) , 
cos(A$ KjK+ i) , 
$ K - , 
cos 9^ , 



(20) 

(21) 
(22) 
(23) 
(24) 
(25) 
(26) 



where the azimuthal scattering angle is A<& K K+ i (equivalent to 4> in the scattering plane representation). The modified 
scattering amplitudes, which depend only on the cosine of the scattering angle /z, due to the spherical symmetry of 
the scatterer, 



X2(p) 



S\\(n) - fJ.S ± (n) 

l-/i 2 

cos 6 = (l,l)(r,r) - (l,r)(r,l) 

\J( l - - ^l+i) cos(A$ K:K+ i) + n K fi K+1 



(27) 
(28) 

(29) 



with the scattering amplitudes given by Eqs.|3and^ Be- 
cause of the basis transformations, the Jones matrix in 
the Chandrasekhar-Sekera representation Eq. H2U|) can- 
not simply be written as product of a scattering by a rota- 
tion matrix as in the scattering plane representation (Eq. 
0). Nevertheless, writing k K — z, K = x, and <3? K = y, 
then <& K = 0, and 9 K = tt/2 for the incident direction and 
$ K+ i = <p, and 6 K+ i = 0, then A$ KjK+1 = -(f), /i K = 1, 



and /i K +i = cos 6, the scattering plane representation can 
be obtained as a particular case. 

Although the Chandrasekhar-Sekera representation is 
more complex than the scattering plane representation, it 
is suitable for numerical calculation in the multiple scat- 
tering context since it already includes the basis trans- 
formation in itself. In this representation, given k K and 
k K+ i in the laboratory frame, the electric field along the 
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scattering sequence can be promptly obtained, even in 
the case where the reverse sequences must be calculated 
such as in the weak localization of eletromagnetic waves. 



are given by Eqs. (|12|) and (|13|) . For small particle scat- 
tering, = ia(ka) 3 /2, which leads to X\ = ia(ka) 3 
and X 2 = 0. 



IV. NON-COMMUTABILITY OF LIMITING 
CASES 

To calculate the scattering of light by polydisper- 
sions, particle size may range from the small particle to 
large sphere limit. Numerical problems with Eqs. I|27() 
and (|28|) may occur since for/backward scattering events 
(/i — * ±1) have non- vanishing probability. 

Back/forward scattering are important in practice. For 
instance, the phase functions are very forwardly peaked 
(/j — >• 1) in the Mie regime. In lidar calculation, the 
interest is in the backward scattering events (/i — ► —1, 
in general these scattering events are forced to occur by 
Metropolis variance-reduction methods in a Monte Carlo 
scheme) . 

Since these difficulties have never been reported in the 
literature, it is explicitly developed here, where we con- 
sider the particles to be non-magnetic to = [i s j [i m = 1. 

A simple way to obtain the limiting values, which 
works properly for the optical theorem, is to rewrite 
Eqs. E7J and JSJJ as: X 1 (fi) = [S±(ji) - fJ,S\\ (/i)]/[(l + 
/i)(l-/x)] and A 2 (/i) = [5||(/i)-M^±(M)]/[(l + M)(l-M)]- 
In the limiting cases of fj, — > ±1, <Sj_(l) = ^(l) = 
S(l), S±(-l) = -S{-1) and = S{-1), with 

S'(±l) given by Eq. l|l(jH an d one can write: Xi(fi — > 
±1) = ±5(±l)/2 and X 2 {p -► ±1) = S*(±l)/2. In 
short, Xi(±\) = ±X2(±1). Now let us consider the 
limit ka — » (small particle limit), the leading term is 
a\ = —ia(ka) 3 . The scattering amplitudes are: Sj_ = 
— ia(ka) 3 (independent of /j) and S\\(fi) = S±(a with 
a = (to 2 -1)/(to 2 +2) 0,1113 so that: Xi((i) = -ia(ka) 3 
and X 2 (fJ.) = 0. This implies Xx(l) = S± = S\\ and 
X 2 (l) = leading to Xi(±l) ^ ±A 2 (±1). Formally, X x 
tends to zero, as ka tends to zero, while X 2 vanishes. 
Taking the limit fj, — ► ±1 and then ka — > the result is 
different than taking ka — » and then \i — ► ±1. These 
limits do not commute, which is a non-physical result. 
Furthermore, the order the limits are taken affects nu- 
merical results. 

Once this problem has been noticed, the solution is ob- 
vious: expand the scattering amplitudes to higher-order 
terms as function of cos 9. Because of the presence of a 
1 — /i 2 = sin 2 9 in the denominator of Xx (//) and X 2 (/i) 
(Eqs. (|2"7Jl and 1 (25) 1 1. the consideration up to fourth- 
order terms in 9 is mandatory here, in contrast to the 
use of the optical theorem. 

For fj, — ► ±1, using Eqs. and and and 

1(15) one obtains: -> ±1) = ±5(±l)/2 + A' 1 ', 

X 2 (v -» ±1) = ^±ii - A(±), with AC±) - ±P(±) - 



1)(2^ ±} 



L<±> = (l/2)E^=i(2n 
6„)/[n(n 4- 1)], which corrects the problem of non- 

commutability of limits. The numbers and 



V. THE NEW SET OF BASIS FUNCTIONS 

An effective manner to solve the non-commutability of 
the limiting cases is to use the definitions of the scattering 
amplitudes (Eqs. © and (0J) and rewrite Xi(/j) and 
X 2 (ii) in Eqs. and as: 



X 2 (p) = 

xl 2) M = 



Y> 2n+l 
2— 1 n(n + 1) 

n— 1 x ' 

A 2n+l 



^-f n(rc+l) 

n— 1 v 7 
7T n (/x) ~ ^ T «(M) 

1-/1 2 
Tn(/x) - fmnJji} 

1-/1 2 



a„xl 2) (M)+rfW) 

= 7r„( M ) + M<M (32) 
= -<(a0 , (33) 



where Eqs. J5J) and © have been used. 

Instead of calculating the scattering amplitudes S± (/i) 
and S\\(n), Xi(fj,) and X 2 (fi) can be directly calculated 

from the basis functions Xn\p) and Xn\^)- The con- 
sideration of Eqs. H32|) and l|33(l implies to an important 
achievement, since the term sin 2 9 that appeared in the 
denominator of Eqs. (|27|l and (|28|l cancels out when one 
writes t„(/i) = /«7r n (//) — (1 — y?)-K' n {y) (Eq. JO)). In this 
way, no asymptotic behavior for Xi(n) and X 2 (fi) must 
be explicitally calculated numerically, avoiding numeri- 
cal tests, which increase computer time. This procedure 
achieves directly the limiting cases, independently of the 
order these limits are taken. 

The functions 7r„(/i) and Tt' n (fi) can be obtained nu- 
merically by stable recursion relations 0, Hfj 



t n -i(p) = (2n- l)7r n _i(/i) 

Tn(M) = ^n-l(^) + <-2(M) 

/ \ M^i-i(m) - nn n _ 2 (ij,) 

M/*) = ; 

n — 1 



(34) 
(35) 

(36) 



which are initialized with 7ro(/x) = 0, 7Ti (/x) = 1, ir 2 (p) = 
3/it, 7Ti(/i) = and 7r 2 (/i) = 3. 

A further improvement can be achieved in numerical 
calculations. Instead of using Xi(/i) and A 2 (/i), one 
writes: 



1 



2n± 1 



2 ^ n(„ + 1) 



(a n ±&„)xW( M ) (37) 



(38) 



which is an adaptation of the Wiscombe's method [39| . 
For small particle scattering A±(/i) = — 3«a(/ca) 3 /4, this 
leads to the correct limits X\ = —ta(ka) 3 and X 2 = 0. 
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In brief the new algorithm works as follows. Given 
the size parameter fca, the relative refraction index m 
and fh (for magnetic scattering), one calculates the Mie 
coefficients a n and b n and store the terms (2n + l)(a n ± 
b n )/[2n(n+l)] (Eq. (|3"?|) ) in arrays |3S| . Given the input 
k K and output fc K +i directions, proceed as follows: 

1. obtain the scattering angle as given by Eq. i|29[l : 

2. calculate the "Legendre" polynomials by the recur- 
sion relationships of Eqs. (|3"4"1 ) l|3~5l) and 

3. the basis functions xt (m) are obtained by Eq. (jSSJ; 

4. A±(^), Xi(p) and X 2 (^) are given by 
Eqs. (??), P7|) and $5T\; 

5. the azimuthal angle is calculated with Eq. I|25[) : 

6. the director cosines are obtained by Eq. I|26|) ; 

7. use Eqs. (Ell), I®, © and (J23J to obtain (Z,Z), 
(Z,r), (r, ?) and (r, r), respectively; 

8. the results of the above items permit the calcula- 
tions of the Jones matrix with Eq. (|20(l . 

9. and finally, at the distance R K +\ from the center of 
the sphere along k K +x, the scattered electric field is 
given by Eq. l|I5fr. 

The calculations are performed on the field level. To 
obtain the Stokes intensities (I, Q, U and V) use the 
coherence matrix E^ +1 E K+ i ^3% . 

For a fixed particle size (monodispersions), the func- 
tions Xn(fJ-) are calculated each time the angles of 
source and detection are modified, or when a new scat- 
tering event occurs in a Monte Carlo sequence. In this 
case, it is suitable to tabulate X±(fi) (Eq. ET7)l . in un- 
even partitions. Higher resolution are required for |/x| 
close to unity |3^|, roughly 9 < ka in the forward case 
and 7r — 8 < (fca) 4//3 in the backwards case. To consider 
polydispersions, the size-distribution function must be 



taken into account. Some of these size-distribution func- 
tions are found in Ref. 47] and an implementation for 
large size parameters is found in Ref. [48j. In this case, 
values of (\X+(^), (|X_(/i)| 2 }, Re[(X + ( M )X* ( M ))] and 
Im[(X + (/i)Xl (/i))] should be tabulated, where (. . .) 
refers to the size-distribution average. 

VI. CONCLUSION 

The Chandrasekhar-Sekera representation is more ap- 
propiate than the scattering plane representation to be 
considered when multiple M ie scattering of light is taken 
into account. Nevertheless care must be taken when this 
representation is used because of the non-commutability 
of the limiting cases of: Rayleigh and for/backwards scat- 
tering. It must be corrected considering the higher-order 
terms in the scattering angle in the scattering amplitude 
expansions. A new set of basis functions for the calcu- 
lation of the angular functions in the Mie scattering in 
the Chandrasekhar-Sekera representation has been pre- 
sented. An effective algorithm has been implemented us- 
ing new basis functions, which have been introduced. 

Furthermore, the use of a generalization (to associated 

Legendre polynomials) of XrP and Xn^ can be used for 
non-spherical scattering, such as in the extended bound- 
ary condition presented in Ref. |49| . This generalization, 
among the generalization for the spherical Henkel func- 
tions, permits the writing of the scattering amplitudes 
of non-spherical particles in terms of matrixes, which are 
suitable for numerical calculations. These calculations 
will be reported in future. 
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